home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / frasrc19.zip / BIGINIT.C < prev    next >
C/C++ Source or Header  |  1995-03-04  |  23KB  |  578 lines
  1. /* biginit.c - C routines for bignumbers */
  2.  
  3. /*
  4. Wesley Loewer's Big Numbers.        (C) 1994, Wesley B. Loewer
  5. */
  6.  
  7. #include <stdlib.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include <malloc.h>
  11. #include <math.h>
  12. #include "fractint.h"
  13. #include "fractype.h"
  14. #include "prototyp.h"
  15. #include "biginit.h"
  16.  
  17. /* appears to me that avoiding the start of extraseg is unnecessary. If
  18.    correct, later we can eliminate ENDVID here. */
  19. #ifdef ENDVID
  20. #undef ENDVID
  21. #endif
  22. #define ENDVID 0
  23.  
  24. /* globals */
  25. int bnstep, bnlength, intlength, rlength, padding, shiftfactor, decimals;
  26. int bflength, rbflength, bfshiftfactor, bfdecimals;
  27.  
  28. /* used internally by bignum.c routines */
  29. static bn_t bnroot=NULL;
  30. static bn_t stack_ptr; /* memory allocator base after global variables */
  31. bn_t bntmp1, bntmp2, bntmp3, bntmp4, bntmp5, bntmp6; /* rlength  */
  32. bn_t bntmpcpy1, bntmpcpy2;                           /* bnlength */
  33.  
  34. /* used by other routines */
  35. bn_t bnxmin, bnxmax, bnymin, bnymax, bnx3rd, bny3rd;        /* bnlength */
  36. bn_t bnxdel, bnydel, bnxdel2, bnydel2, bnclosenuff;         /* bnlength */
  37. bn_t bntmpsqrx, bntmpsqry, bntmp;                           /* rlength  */
  38. _BNCMPLX bnold, /* bnnew, */ bnparm, bnsaved;               /* bnlength */
  39. _BNCMPLX bnnew;                                              /* rlength */
  40. bn_t bn_pi;                                           /* TAKES NO SPACE */
  41.  
  42. bf_t bftmp1, bftmp2, bftmp3, bftmp4, bftmp5, bftmp6;     /* rbflength+2 */
  43. bf_t bftmpcpy1, bftmpcpy2;                               /* rbflength+2 */
  44. bf_t bfxdel, bfydel, bfxdel2, bfydel2, bfclosenuff;      /* rbflength+2 */
  45. bf_t bftmpsqrx, bftmpsqry;                               /* rbflength+2 */
  46. _BFCMPLX /* bfold,  bfnew, */ bfparm, bfsaved;            /* bflength+2 */
  47. _BFCMPLX bfold,  bfnew;                                  /* rbflength+2 */
  48. bf_t bf_pi;                                           /* TAKES NO SPACE */
  49. bf_t big_pi;                                              /* bflength+2 */
  50.  
  51. /* for testing only */
  52.  
  53. /* used by other routines */
  54. bf_t bfxmin, bfxmax, bfymin, bfymax, bfx3rd, bfy3rd;      /* bflength+2 */
  55. bf_t bfsxmin, bfsxmax, bfsymin, bfsymax, bfsx3rd, bfsy3rd;/* bflength+2 */
  56. bf_t bfparms[10];                                    /* (bflength+2)*10 */
  57. bf_t bftmp;
  58.  
  59. #define LOG10_256 2.4082399653118
  60. #define LOG_256   5.5451774444795
  61.  
  62. static int save_bf_vars(void);
  63. static int restore_bf_vars(void);
  64.  
  65. #if CALCULATING_BIG_PI  /* normally used */
  66. static void generate_big_pi(void);
  67. #endif
  68.  
  69. /*********************************************************************/
  70. /* given bnlength, calc_lengths will calculate all the other lengths */
  71. void calc_lengths(void)
  72.     {
  73. #if 0
  74. #ifdef USE_BIGNUM_C_CODE
  75.     bnstep = 2;
  76. #else /* use 80x86 asm code */
  77.     if (cpu >= 386)
  78.         bnstep = 4;
  79.     else /* cpu <= 286 */
  80.         bnstep = 2;
  81. #endif
  82. #else
  83.         bnstep = 4;  /* use 4 in all cases */
  84. #endif
  85.  
  86.     if (bnlength % bnstep != 0)
  87.         bnlength = (bnlength / bnstep + 1) * bnstep;
  88.     if (bnlength == bnstep)
  89.         padding = bnlength;
  90.     else
  91.         padding = 2*bnstep;
  92.     rlength = bnlength + padding;
  93.  
  94.     /* This shiftfactor assumes non-full multiplications will be performed.*/
  95.     /* Change to bnlength-intlength for full multiplications.              */
  96.     shiftfactor = padding - intlength;
  97.  
  98.     bflength = bnlength+bnstep; /* one extra step for added precision */
  99.     rbflength = bflength + padding;
  100.     bfdecimals = (int)((bflength-2)*LOG10_256);
  101.     bfshiftfactor = padding - 2;
  102.     }
  103.  
  104. /************************************************************************/
  105. /* intended only to be called from init_bf_dec() or init_bf_length().   */
  106. /* initialize bignumber global variables                                */
  107.  
  108. long maxptr = 0;
  109. long startstack = 0;
  110. long maxstack = 0;
  111. int bf_save_len = 0;
  112.  
  113. /* ??? for some strange reason, msc 7.0 hangs here without this pragma. ??? */
  114. #pragma optimize( "", off )
  115. static void init_bf_2(void)
  116.     {
  117.     int i;
  118.     long ptr;
  119.     save_bf_vars(); /* copy corners values for conversion */
  120.  
  121.     calc_lengths();
  122.  
  123.     /* allocate all the memory at once within the same segment (DOS) */
  124.     bnroot = (bf_t)MK_FP(extraseg,ENDVID);  /* ENDVID is to avoid videotable */
  125.  
  126.     /* at present time one call would suffice, but this logic allows
  127.        mulytiple kinds of alternate math eg long double */ 
  128.     if((i = find_alternate_math(fractype, BIGNUM)) > -1)
  129.        bf_math = alternatemath[i].math;
  130.     else if((i = find_alternate_math(fractype, BIGFLT)) > -1)
  131.        bf_math = alternatemath[i].math;
  132.     else
  133.        bf_math = 1; /* maybe called from cmdfiles.c and fractype not set */
  134.  
  135.     floatflag=1;
  136.  
  137.     /* Now split up the memory among the pointers */
  138.     /* internal pointers */
  139.     ptr        = 0;
  140.     bntmp1     = bnroot+ptr; ptr += rlength;
  141.     bntmp2     = bnroot+ptr; ptr += rlength;  
  142.     bntmp3     = bnroot+ptr; ptr += rlength;
  143.     bntmp4     = bnroot+ptr; ptr += rlength;
  144.     bntmp5     = bnroot+ptr; ptr += rlength;
  145.     bntmp6     = bnroot+ptr; ptr += rlength;
  146.  
  147.     bftmp1     = bnroot+ptr; ptr += rbflength+2;
  148.     bftmp2     = bnroot+ptr; ptr += rbflength+2;
  149.     bftmp3     = bnroot+ptr; ptr += rbflength+2;
  150.     bftmp4     = bnroot+ptr; ptr += rbflength+2;
  151.     bftmp5     = bnroot+ptr; ptr += rbflength+2;
  152.     bftmp6     = bnroot+ptr; ptr += rbflength+2;
  153.  
  154.     bftmpcpy1  = bnroot+ptr; ptr += (rbflength+2)*2;
  155.     bftmpcpy2  = bnroot+ptr; ptr += (rbflength+2)*2;
  156.  
  157.     bntmpcpy1  = bnroot+ptr; ptr += (rlength*2);
  158.     bntmpcpy2  = bnroot+ptr; ptr += (rlength*2);
  159.  
  160.     if (bf_math == BIGNUM)
  161.     {
  162.     bnxmin     = bnroot+ptr; ptr += bnlength;
  163.     bnxmax     = bnroot+ptr; ptr += bnlength;
  164.     bnymin     = bnroot+ptr; ptr += bnlength;
  165.     bnymax     = bnroot+ptr; ptr += bnlength;
  166.     bnx3rd     = bnroot+ptr; ptr += bnlength;
  167.     bny3rd     = bnroot+ptr; ptr += bnlength;
  168.     bnxdel     = bnroot+ptr; ptr += bnlength;
  169.     bnydel     = bnroot+ptr; ptr += bnlength;
  170.     bnxdel2    = bnroot+ptr; ptr += bnlength;
  171.     bnydel2    = bnroot+ptr; ptr += bnlength;
  172.     bnold.x    = bnroot+ptr; ptr += rlength;
  173.     bnold.y    = bnroot+ptr; ptr += rlength;
  174.     bnnew.x    = bnroot+ptr; ptr += rlength;
  175.     bnnew.y    = bnroot+ptr; ptr += rlength;
  176.     bnsaved.x  = bnroot+ptr; ptr += bnlength;
  177.     bnsaved.y  = bnroot+ptr; ptr += bnlength;
  178.     bnclosenuff= bnroot+ptr; ptr += bnlength;
  179.     bnparm.x   = bnroot+ptr; ptr += bnlength;
  180.     bnparm.y   = bnroot+ptr; ptr += bnlength;
  181.     bntmpsqrx  = bnroot+ptr; ptr += rlength;
  182.     bntmpsqry  = bnroot+ptr; ptr += rlength;
  183.     bntmp      = bnroot+ptr; ptr += rlength;
  184.     }
  185.     if (bf_math == BIGFLT)
  186.     {
  187.     bfxdel     = bnroot+ptr; ptr += bflength+2;
  188.     bfydel     = bnroot+ptr; ptr += bflength+2;
  189.     bfxdel2    = bnroot+ptr; ptr += bflength+2;
  190.     bfydel2    = bnroot+ptr; ptr += bflength+2;
  191.     bfold.x    = bnroot+ptr; ptr += rbflength+2;
  192.     bfold.y    = bnroot+ptr; ptr += rbflength+2;
  193.     bfnew.x    = bnroot+ptr; ptr += rbflength+2;
  194.     bfnew.y    = bnroot+ptr; ptr += rbflength+2;
  195.     bfsaved.x  = bnroot+ptr; ptr += bflength+2;
  196.     bfsaved.y  = bnroot+ptr; ptr += bflength+2;
  197.     bfclosenuff= bnroot+ptr; ptr += bflength+2;
  198.     bfparm.x   = bnroot+ptr; ptr += bflength+2;
  199.     bfparm.y   = bnroot+ptr; ptr += bflength+2;
  200.     bftmpsqrx  = bnroot+ptr; ptr += rbflength+2;
  201.     bftmpsqry  = bnroot+ptr; ptr += rbflength+2;
  202.     big_pi     = bnroot+ptr; ptr += bflength+2;
  203.     bftmp      = bnroot+ptr; ptr += rbflength+2;
  204.     }
  205.     stack_ptr  = bnroot + ptr;    
  206.     startstack = ptr;
  207.  
  208.     /* max stack offset from bnroot */
  209.     maxstack = (long)0x10000l-(bflength+2)*22-ENDVID;
  210.  
  211.     /* sanity check */ 
  212.     /* leave room for NUMVARS variables allocated from stack */
  213.     /* also leave room for the safe area at top of segment */
  214.     if(ptr + NUMVARS*(bflength+2) > (long)0x10000l -(bflength+2)*22-ENDVID)
  215.        {
  216.        char msg[80];
  217.        char nmsg[80];
  218.        static FCODE fmsg[] = {"Requested precision of %d too high, aborting"};
  219.        far_strcpy(nmsg,fmsg);
  220.        sprintf(msg,nmsg,decimals);
  221.        stopmsg(0,msg);
  222.        goodbye();
  223.        }
  224.  
  225.     /* room for 6 corners + 6 save corners + 10 params at top of extraseg */
  226.     /* this area is safe - use for variables that are used outside fractal*/
  227.     /* generation - e.g. zoom box variables */
  228.     ptr  = ((long)0x10000l-(bflength+2)*22-ENDVID);
  229.     bfxmin     = bnroot+ptr; ptr += bflength+2;
  230.     bfxmax     = bnroot+ptr; ptr += bflength+2;
  231.     bfymin     = bnroot+ptr; ptr += bflength+2;
  232.     bfymax     = bnroot+ptr; ptr += bflength+2;
  233.     bfx3rd     = bnroot+ptr; ptr += bflength+2;
  234.     bfy3rd     = bnroot+ptr; ptr += bflength+2;
  235.     for(i=0;i<10;i++) 
  236.        {
  237.        bfparms[i]  = bnroot+ptr; ptr += bflength+2;
  238.        }    
  239.     bfsxmin    = bnroot+ptr; ptr += bflength+2;
  240.     bfsxmax    = bnroot+ptr; ptr += bflength+2;
  241.     bfsymin    = bnroot+ptr; ptr += bflength+2;
  242.     bfsymax    = bnroot+ptr; ptr += bflength+2;
  243.     bfsx3rd    = bnroot+ptr; ptr += bflength+2;
  244.     bfsy3rd    = bnroot+ptr; ptr += bflength+2;
  245.     /* end safe vars */     
  246.  
  247.     /* good citizens initialize variables */
  248.     if(bf_save_len)  /* leave save area */
  249.        far_memset(bnroot+(bf_save_len+2)*22,0,(unsigned)(startstack-(bf_save_len+2)*22));
  250.     else /* first time through - nothing saved */
  251.        {
  252.        /* high variables */
  253.        far_memset(bnroot+((long)0x10000l-(bflength+2)*22-ENDVID),0,(bflength+2)*22);
  254.        /* low variables */
  255.        far_memset(bnroot,0,(unsigned)startstack);
  256.        } 
  257.  
  258.     restore_bf_vars();
  259.  
  260. #if CALCULATING_BIG_PI  /* not normally used */
  261.     generate_big_pi();
  262. #endif
  263.  
  264.     /* Initialize the value of pi.  Needed for trig functions. */
  265.     /* init_big_pi(); */
  266. /* call to init_big_pi() has been moved to fractal setup routine */
  267. /* so as to use only when necessary. */
  268.  
  269.     }
  270.  
  271.  
  272. /**********************************************************/
  273. /* save current corners and parameters to start of bnroot */
  274. /* to preserve values across calls to init_bf()           */
  275. static int save_bf_vars(void)
  276.    {
  277.    int ret;
  278.    unsigned int mem;
  279.    if(bnroot != NULL)
  280.       {
  281.       mem = (bflength+2)*22;  /* 6 corners + 6 save corners + 10 params */
  282.       bf_save_len = bflength;
  283.       far_memcpy(bnroot,bfxmin,mem);
  284.       /* scrub old high area */
  285.       far_memset(bfxmin,0,mem);
  286.       ret = 0;
  287.       }
  288.    else 
  289.       {
  290.       bf_save_len = 0;   
  291.       ret = -1;
  292.       }
  293.    return(ret);   
  294.    }
  295.  
  296. /************************************************************************/
  297. /* copy current corners and parameters from save location               */
  298. static int restore_bf_vars(void)
  299.    {
  300.    bf_t ptr;
  301.    int i;
  302.    if(bf_save_len == 0)
  303.       return(-1);
  304.    ptr  = bnroot;
  305.    convert_bf(bfxmin,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  306.    convert_bf(bfxmax,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  307.    convert_bf(bfymin,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  308.    convert_bf(bfymax,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  309.    convert_bf(bfx3rd,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  310.    convert_bf(bfy3rd,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  311.    for(i=0;i<10;i++)
  312.       {
  313.       convert_bf(bfparms[i],ptr,bflength,bf_save_len);
  314.       ptr += bf_save_len+2;
  315.       }     
  316.    convert_bf(bfsxmin,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  317.    convert_bf(bfsxmax,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  318.    convert_bf(bfsymin,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  319.    convert_bf(bfsymax,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  320.    convert_bf(bfsx3rd,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  321.    convert_bf(bfsy3rd,ptr,bflength,bf_save_len); ptr += bf_save_len+2;
  322.  
  323.    /* scrub save area */
  324.    far_memset(bnroot,0,(bf_save_len+2)*22);
  325.    return(0);
  326.    }
  327.  
  328. /*******************************************/
  329. /* free corners and parameters save memory */
  330. void free_bf_vars()
  331.    {
  332.    bf_save_len = bf_math = 0;
  333.    bnstep=bnlength=intlength=rlength=padding=shiftfactor=decimals=0;
  334.    bflength=rbflength=bfshiftfactor=bfdecimals=0;
  335.    }
  336.       
  337. #pragma optimize( "", on )
  338.  
  339.  
  340. /************************************************************************/
  341. /* Memory allocator routines start here.                                */
  342. /************************************************************************/
  343. /* Allocates a bn_t variable on stack                                   */
  344. bn_t alloc_stack(size_t size)
  345.    {
  346.    long stack_addr;
  347.    if(bf_math == 0)
  348.       {
  349.       static FCODE msg[] = {"alloc_stack called with bf_math==0"};
  350.       stopmsg(0,msg);
  351.       return(0);
  352.       }
  353.    stack_addr = (long)(stack_ptr-bnroot)+size+ENDVID;
  354.    
  355.    if((stack_addr) > ((long)0x10000l-(bflength+2)*22))
  356.       {
  357.       static FCODE msg[] = {"Aborting, Out of Bignum Stack Space"};
  358.       stopmsg(0,msg);
  359.       goodbye();     
  360.       }
  361.    /* keep track of max ptr */
  362.    if(stack_addr > maxptr) 
  363.       maxptr = stack_addr;
  364.    stack_ptr += size;   /* increment stack pointer */
  365.    return(stack_ptr - size);
  366.    }
  367.    
  368. /************************************************************************/
  369. /* Returns stack pointer offset so it can be saved.                            */
  370. int save_stack(void)
  371.    {
  372.    return(stack_ptr - bnroot);
  373.    }
  374.  
  375. /************************************************************************/
  376. /* Restores stack pointer, effectively freeing local variables          */
  377. /*    allocated since save_stack()                                   */
  378. void restore_stack(int old_offset)
  379.    {
  380.    stack_ptr  = bnroot+old_offset;
  381.    }
  382.  
  383. /************************************************************************/
  384. /* Memory allocator routines end here.                                  */
  385. /************************************************************************/
  386.  
  387. /************************************************************************/
  388. /* initialize bignumber global variables                                */
  389. /*   dec = decimal places after decimal point                           */
  390. /*   intl = bytes for integer part (1, 2, or 4)                         */
  391.  
  392. void init_bf_dec(int dec)
  393.     {
  394.     if(bfdigits)
  395.        decimals=bfdigits;   /* blindly force */
  396.     else   
  397.        decimals = dec;
  398.     if (fractype == FPMANDELZPOWER || fractype == FPJULIAZPOWER)
  399.        intlength = 2;
  400.     /* the bailout tests need greater dynamic range */
  401.     else if(bailoutest == Real || bailoutest == Imag || bailoutest == And)
  402.        intlength = 2;
  403.     else
  404.        intlength = 1;   
  405.     /* conservative estimate */
  406.     bnlength = intlength + (int)(decimals/LOG10_256) + 1; /* round up */
  407.     init_bf_2();
  408.     }
  409.  
  410. /************************************************************************/
  411. /* initialize bignumber global variables                                */
  412. /*   bnl = bignumber length                                             */
  413. /*   intl = bytes for integer part (1, 2, or 4)                         */
  414. void init_bf_length(int bnl)
  415.     {
  416.     bnlength = bnl;
  417.  
  418.     if (fractype == FPMANDELZPOWER || fractype == FPJULIAZPOWER)
  419.        intlength = 2;
  420.     /* the bailout tests need greater dynamic range */
  421.     else if(bailoutest == Real || bailoutest == Imag || bailoutest == And)
  422.        intlength = 2;
  423.     else
  424.        intlength = 1;   
  425.     /* conservative estimate */
  426.     decimals = (int)((bnlength-intlength)*LOG10_256);
  427.     init_bf_2();
  428.     }
  429.  
  430.  
  431. void init_big_pi(void)
  432.     {
  433.     /* What, don't you recognize the first 700 digits of pi, */
  434.     /* in base 256, in reverse order?                        */
  435.     int length, pi_offset;
  436.     static BFCODE pi_table[] = {
  437.          0x44, 0xD5, 0xDB, 0x69, 0x17, 0xDF, 0x2E, 0x56, 0x87, 0x1A,
  438.          0xA0, 0x8C, 0x6F, 0xCA, 0xBB, 0x57, 0x5C, 0x9E, 0x82, 0xDF,
  439.          0x00, 0x3E, 0x48, 0x7B, 0x31, 0x53, 0x60, 0x87, 0x23, 0xFD,
  440.          0xFA, 0xB5, 0x3D, 0x32, 0xAB, 0x52, 0x05, 0xAD, 0xC8, 0x1E,
  441.          0x50, 0x2F, 0x15, 0x6B, 0x61, 0xFD, 0xDF, 0x16, 0x75, 0x3C,
  442.          0xF8, 0x22, 0x32, 0xDB, 0xF8, 0xE9, 0xA5, 0x8E, 0xCC, 0xA3,
  443.          0x1F, 0xFB, 0xFE, 0x25, 0x9F, 0x67, 0x79, 0x72, 0x2C, 0x40,
  444.          0xC6, 0x00, 0xA1, 0xD6, 0x0A, 0x32, 0x60, 0x1A, 0xBD, 0xC0,
  445.          0x79, 0x55, 0xDB, 0xFB, 0xD3, 0xB9, 0x39, 0x5F, 0x0B, 0xD2,
  446.          0x0F, 0x74, 0xC8, 0x45, 0x57, 0xA8, 0xCB, 0xC0, 0xB3, 0x4B,
  447.          0x2E, 0x19, 0x07, 0x28, 0x0F, 0x66, 0xFD, 0x4A, 0x33, 0xDE,
  448.          0x04, 0xD0, 0xE3, 0xBE, 0x09, 0xBD, 0x5E, 0xAF, 0x44, 0x45,
  449.          0x81, 0xCC, 0x2C, 0x95, 0x30, 0x9B, 0x1F, 0x51, 0xFC, 0x6D,
  450.          0x6F, 0xEC, 0x52, 0x3B, 0xEB, 0xB2, 0x39, 0x13, 0xB5, 0x53,
  451.          0x6C, 0x3E, 0xAF, 0x6F, 0xFB, 0x68, 0x63, 0x24, 0x6A, 0x19,
  452.          0xC2, 0x9E, 0x5C, 0x5E, 0xC4, 0x60, 0x9F, 0x40, 0xB6, 0x4F,
  453.          0xA9, 0xC1, 0xBA, 0x06, 0xC0, 0x04, 0xBD, 0xE0, 0x6C, 0x97,
  454.          0x3B, 0x4C, 0x79, 0xB6, 0x1A, 0x50, 0xFE, 0xE3, 0xF7, 0xDE,
  455.          0xE8, 0xF6, 0xD8, 0x79, 0xD4, 0x25, 0x7B, 0x1B, 0x99, 0x80,
  456.          0xC9, 0x72, 0x53, 0x07, 0x9B, 0xC0, 0xF1, 0x49, 0xD3, 0xEA,
  457.          0x0F, 0xDB, 0x48, 0x12, 0x0A, 0xD0, 0x24, 0xD7, 0xD0, 0x37,
  458.          0x3D, 0x02, 0x9B, 0x42, 0x72, 0xDF, 0xFE, 0x1B, 0x06, 0x77,
  459.          0x3F, 0x36, 0x62, 0xAA, 0xD3, 0x4E, 0xA6, 0x6A, 0xC1, 0x56,
  460.          0x9F, 0x44, 0x1A, 0x40, 0x73, 0x20, 0xC1, 0x85, 0xD8, 0x75,
  461.          0x6F, 0xE0, 0xBE, 0x5E, 0x8B, 0x3B, 0xC3, 0xA5, 0x84, 0x7D,
  462.          0xB4, 0x9F, 0x6F, 0x45, 0x19, 0x86, 0xEE, 0x8C, 0x88, 0x0E,
  463.          0x43, 0x82, 0x3E, 0x59, 0xCA, 0x66, 0x76, 0x01, 0xAF, 0x39,
  464.          0x1D, 0x65, 0xF1, 0xA1, 0x98, 0x2A, 0xFB, 0x7E, 0x50, 0xF0,
  465.          0x3B, 0xBA, 0xE4, 0x3B, 0x7A, 0x13, 0x6C, 0x0B, 0xEF, 0x6E,
  466.          0xA3, 0x33, 0x51, 0xAB, 0x28, 0xA7, 0x0F, 0x96, 0x68, 0x2F,
  467.          0x54, 0xD8, 0xD2, 0xA0, 0x51, 0x6A, 0xF0, 0x88, 0xD3, 0xAB,
  468.          0x61, 0x9C, 0x0C, 0x67, 0x9A, 0x6C, 0xE9, 0xF6, 0x42, 0x68,
  469.          0xC6, 0x21, 0x5E, 0x9B, 0x1F, 0x9E, 0x4A, 0xF0, 0xC8, 0x69,
  470.          0x04, 0x20, 0x84, 0xA4, 0x82, 0x44, 0x0B, 0x2E, 0x39, 0x42,
  471.          0xF4, 0x83, 0xF3, 0x6F, 0x6D, 0x0F, 0xC5, 0xAC, 0x96, 0xD3,
  472.          0x81, 0x3E, 0x89, 0x23, 0x88, 0x1B, 0x65, 0xEB, 0x02, 0x23,
  473.          0x26, 0xDC, 0xB1, 0x75, 0x85, 0xE9, 0x5D, 0x5D, 0x84, 0xEF,
  474.          0x32, 0x80, 0xEC, 0x5D, 0x60, 0xAC, 0x7C, 0x48, 0x91, 0xA9,
  475.          0x21, 0xFB, 0xCC, 0x09, 0xD8, 0x61, 0x93, 0x21, 0x28, 0x66,
  476.          0x1B, 0xE8, 0xBF, 0xC4, 0xAF, 0xB9, 0x4B, 0x6B, 0x98, 0x48,
  477.          0x8F, 0x3B, 0x77, 0x86, 0x95, 0x28, 0x81, 0x53, 0x32, 0x7A,
  478.          0x5C, 0xCF, 0x24, 0x6C, 0x33, 0xBA, 0xD6, 0xAF, 0x1E, 0x93,
  479.          0x87, 0x9B, 0x16, 0x3E, 0x5C, 0xCE, 0xF6, 0x31, 0x18, 0x74,
  480.          0x5D, 0xC5, 0xA9, 0x2B, 0x2A, 0xBC, 0x6F, 0x63, 0x11, 0x14,
  481.          0xEE, 0xB3, 0x93, 0xE9, 0x72, 0x7C, 0xAF, 0x86, 0x54, 0xA1,
  482.          0xCE, 0xE8, 0x41, 0x11, 0x34, 0x5C, 0xCC, 0xB4, 0xB6, 0x10,
  483.          0xAB, 0x2A, 0x6A, 0x39, 0xCA, 0x55, 0x40, 0x14, 0xE8, 0x63,
  484.          0x62, 0x98, 0x48, 0x57, 0x94, 0xAB, 0x55, 0xAA, 0xF3, 0x25,
  485.          0x55, 0xE6, 0x60, 0x5C, 0x60, 0x55, 0xDA, 0x2F, 0xAF, 0x78,
  486.          0x27, 0x4B, 0x31, 0xBD, 0xC1, 0x77, 0x15, 0xD7, 0x3E, 0x8A,
  487.          0x1E, 0xB0, 0x8B, 0x0E, 0x9E, 0x6C, 0x0E, 0x18, 0x3A, 0x60,
  488.          0xB0, 0xDC, 0x79, 0x8E, 0xEF, 0x38, 0xDB, 0xB8, 0x18, 0x79,
  489.          0x41, 0xCA, 0xF0, 0x85, 0x60, 0x28, 0x23, 0xB0, 0xD1, 0xC5,
  490.          0x13, 0x60, 0xF2, 0x2A, 0x39, 0xD5, 0x30, 0x9C, 0xB5, 0x59,
  491.          0x5A, 0xC2, 0x1D, 0xA4, 0x54, 0x7B, 0xEE, 0x4A, 0x15, 0x82,
  492.          0x58, 0xCD, 0x8B, 0x71, 0x58, 0xB6, 0x8E, 0x72, 0x8F, 0x74,
  493.          0x95, 0x0D, 0x7E, 0x3D, 0x93, 0xF4, 0xA3, 0xFE, 0x58, 0xA4,
  494.          0x69, 0x4E, 0x57, 0x71, 0xD8, 0x20, 0x69, 0x63, 0x16, 0xFC,
  495.          0x8E, 0x85, 0xE2, 0xF2, 0x01, 0x08, 0xF7, 0x6C, 0x91, 0xB3,
  496.          0x47, 0x99, 0xA1, 0x24, 0x99, 0x7F, 0x2C, 0xF1, 0x45, 0x90,
  497.          0x7C, 0xBA, 0x96, 0x7E, 0x26, 0x6A, 0xED, 0xAF, 0xE1, 0xB8,
  498.          0xB7, 0xDF, 0x1A, 0xD0, 0xDB, 0x72, 0xFD, 0x2F, 0xAC, 0xB5,
  499.          0xDF, 0x98, 0xA6, 0x0B, 0x31, 0xD1, 0x1B, 0xFB, 0x79, 0x89,
  500.          0xD9, 0xD5, 0x16, 0x92, 0x17, 0x09, 0x47, 0xB5, 0xB5, 0xD5,
  501.          0x84, 0x3F, 0xDD, 0x50, 0x7C, 0xC9, 0xB7, 0x29, 0xAC, 0xC0,
  502.          0x6C, 0x0C, 0xE9, 0x34, 0xCF, 0x66, 0x54, 0xBE, 0x77, 0x13,
  503.          0xD0, 0x38, 0xE6, 0x21, 0x28, 0x45, 0x89, 0x6C, 0x4E, 0xEC,
  504.          0x98, 0xFA, 0x2E, 0x08, 0xD0, 0x31, 0x9F, 0x29, 0x22, 0x38,
  505.          0x09, 0xA4, 0x44, 0x73, 0x70, 0x03, 0x2E, 0x8A, 0x19, 0x13,
  506.          0xD3, 0x08, 0xA3, 0x85, 0x88, 0x6A, 0x3F, 0x24,
  507.          /* . */  0x03, 0x00, 0x00, 0x00
  508.               /*  <- up to intlength 4 -> */
  509.               /* or bf_t int length of 2 + 2 byte exp  */
  510.          };
  511.  
  512.     length = bflength+2; /* 2 byte exp */
  513.     pi_offset = sizeof pi_table - length;
  514.     _fmemcpy(big_pi, pi_table + pi_offset, length);
  515.  
  516.     /* notice that bf_pi and bn_pi can share the same memory space */
  517.     bf_pi = big_pi;
  518.     bn_pi = big_pi + (bflength-2) - (bnlength-intlength);
  519.     return;
  520.     }
  521.  
  522.  
  523. #if CALCULATING_BIG_PI
  524. /***********************************************************************/
  525. /* This is a sort of recursive process.  To calculate pi/4,            */
  526. /* atan_bf() will be called, which in turn calls sincos_bf(), which in */
  527. /* turn uses pi/4.  The CALCULATING_BIG_PI is set to 1                 */
  528. /* so that it won't get used by sincos_bf().                           */
  529. /* uses bftmp1 - bftmp6 - global temp bigfloats                        */
  530.  
  531. void generate_big_pi(void)
  532.     {
  533.  
  534.     FILE *f;
  535.     int i, c;
  536.     char str[1710];
  537.     char *s;
  538.  
  539.     /* the following code is used to generate the table of pi */
  540.     stopmsg(0,"Pi is about to be calculated.  This may take a while, so be patient.");
  541.     bf_pi = big_pi;
  542.     inttobf(bf_pi,1);   /* set bf_pi to some value greater than pi/4 */
  543.     inttobf(bftmp6,1);
  544.     unsafe_atan_bf(bf_pi,bftmp6);  /* pi/4 = atan(1) */
  545.     double_a_bf(bf_pi);
  546.     double_a_bf(bf_pi);
  547.  
  548.     /* print in base 256 in reverse order */
  549.     f = fopen("pi256.dat", "wt");
  550.     for (i=0; i<bflength; )
  551.        {
  552.        fprintf(f,"0x%.2X, ", bf_pi[i]);
  553.        if (++i % 10 == 0)
  554.           fprintf(f,"\n");
  555.        }
  556.     fprintf(f,"\n");
  557.     fclose(f);
  558.  
  559.     /* print in base 10 in normal order */
  560.     f = fopen("pi10.dat", "wt");
  561.     bftostr_f(str, 0, bf_pi);
  562.     s = str;
  563.     fprintf(f,"%c", *s++);
  564.     fprintf(f,"%c\n", *s++);
  565.     i = 0;
  566.     while (*s != '\0')
  567.        {
  568.        fprintf(f,"%c", *s++);
  569.        if (++i % 50 == 0)
  570.           fprintf(f,"\n");
  571.        }
  572.     fclose(f);
  573.     stopmsg(0,"Pi calculated and stored in PI256.DAT and PI10.DAT");
  574.     goodbye();
  575.  
  576.     }
  577. #endif
  578.